Trajectory Planning and Simulation Study of Redundant Robotic Arm for Upper Limb Rehabilitation Based on Back Propagation Neural Network and Genetic Algorithm

In this study, a Back Propagation (BP) neural network algorithm based on Genetic Algorithm (GA) optimization is proposed to plan and optimize the trajectory of a redundant robotic arm for the upper limb rehabilitation of patients. The feasibility of the trajectory was verified by numerical simulations. First, the collected dataset was used to train the BP neural network optimized by the GA. Subsequently, the critical points designated by the rehabilitation physician for the upper limb rehabilitation were used as interpolation points for cubic B−spline interpolation to plan the motion trajectory. The GA optimized the planned trajectory with the goal of time minimization, and the feasibility of the optimized trajectory was analyzed with MATLAB simulations. The planned trajectory was smooth and continuous. There was no abrupt change in location or speed. Finally, simulations revealed that the optimized trajectory reduced the motion time and increased the motion speed between two adjacent critical points which improved the rehabilitation effect and can be applied to patients with different needs, which has high application value.


Introduction
Multi−Degree−Of−Freedom (MDOF) robotic arms and automation systems have been widely used in the field of welding and the tools used for welding were often installed on the end of the robotic arm, to replace the worker to complete some complex welding process. The welding action of the robotic arm and the corresponding various parameters have a great impact on the welding effect of the robotic arm. Kang et al. [1] proposed a control algorithm for weaving welding of circular trajectories based on the principle of spatial transformation, which can effectively solve the problem of discontinuities at the sharp corners of the weld, and found that the strategy can generate more accurate circular trajectories. Kang et al. [2] optimized the important parameters of the high−frequency welding process through finite element simulation. In addition, an automated chemical synthesis system was developed to perform chemical experiments [3]. However, the application of robotic arms in the biomedical field has only recently emerged. Bodner et al. [4] experimentally evaluated the suitability of the da Vinci surgical robot for thoracic surgery, proving that it was possible to safely perform thoracic surgery.
There are many factors that can cause damage to the upper limb. For example, stroke [5], prolonged smartphone use [6], and sports−related injuries [7]. Our upper limb rehabilitation robotic arm is not only for a specific type of disease, but also for all patients who need upper limb rehabilitation. Stroke is only used as a typical case to be illustrated in this study. Patients who have had a stroke are difficult to cure in a short time, and most of them require continuous care for rehabilitation [8]. This places a high demand on the number of health care workers. However, most countries have limited medical personnel it is necessary to ensure that the location or speed of the joints did not change abruptly, which eliminated the large errors in the robotic arm and the patients were prevented from receiving additional injuries. That is why we need to perform trajectory planning [22]. Tian et al. [23] proposed a GA to search for effective and optimal solutions in the task trajectory. Gasparetto et al. [24] used a fifth−order B−spline curve for motion trajectory planning and verified the effectiveness of the method via simulation. Sometimes the planned trajectory does not meet the intended requirements, which requires optimization of the trajectory to meet the requirements. Ahuactzin et al. [25] used a GA to optimize the planned trajectory. Therefore, in this study, we planned the rehabilitation trajectory of the robotic arm with the cubic B−spline interpolation method and optimized the planned trajectory by the GA.
To ensure the accuracy and smoothness of the rehabilitation process of the upper limb, the eight critical points in the rehabilitation process were collected from rehabilitation physicians to plan the motion trajectory of the end of the 7DOF redundant upper limb rehabilitation of the robotic arm. To improve the accuracy of the inverse kinematic solution of the robotic arm, as well as the trajectory planning, an algorithm of Back Propagation (BP) neural network optimized by GA was used to calculate the inverse kinematic solution of the robotic arm. The results were used for trajectory planning with the cubic B−spline interpolation method and the trajectory was optimized on the basis of GA for time minimization, which increased the speed of motion between two adjacent critical points and improved the effectiveness of rehabilitation.

Back Propagation Neural Network Algorithm Based on Genetic Algorithm Optimization
The eight critical points were selected by the rehabilitator for the patient's rehabilitation process in the same two−dimensional plane. The motion trajectory consisting of these eight critical points did not guarantee the continuity of speed at the end of the rehabilitation robotic arm, so the corresponding trajectory planning was required to ensure its rehabilitation effect.
In this section, the iterative algorithm of the GA optimized BP neural network was used to find the inverse kinematic solution of the robotic arm. Figure 1 shows the sketch of the robotic arm model. The BP neural network algorithm has the advantages of simple structure and easy implementation. However, it is time−consuming training, easy to fall into local optimum, and no ideal regulation to follow in the selection of the number of nodes in the hidden layer, which makes it impossible to determine the appropriate network structure. Therefore, the weights and thresholds of the neural network was iteratively calculating by a GA. The calculated results were then substituted into the neural network, and the collected dataset was used to train the neural network to acquire the optimal neural network. Finally, the above problem was solved using the optimal network for inverse kinematic computation.
We have fabricated and built a redundant robotic arm for upper limb rehabilitation and have performed some test works on it. Figure 2 shows the SolidWorks model and physical object. The following studies are all based on the upper limb rehabilitation of the 7DOF redundant robotic arm. The BP neural network is a multilayer feedforward neural network trained according to a back−transmission error algorithm, and Figure 3 schematically shows the structure of the neural network of the input layer, hidden layer, and output layer. Each layer contains many neurons inside, and the computational unit in the neuron is the activation function. The activation functions commonly used in network computation are the Sigmoid function [26] and Tansig function [27]. During the operation of the network, data are transmitted sequentially with the input layer, hidden layer, and output layer. When the data reach the output layer they are then transmitted back, thus correcting errors in the network. The network calculation process contains two types of parameters: weights and thresholds. The magnitude of the values affects the error of the network. We have fabricated and built a redundant robotic arm for uppe and have performed some test works on it. Figure 2 shows the So physical object. The following studies are all based on the upper limb 7DOF redundant robotic arm. We have fabricated and built a redundant robotic arm for upper limb rehabilitation and have performed some test works on it. Figure 2 shows the SolidWorks model and physical object. The following studies are all based on the upper limb rehabilitation of the 7DOF redundant robotic arm.

Back Propagation Neural Network
The BP neural network is a multilayer feedforward neural network trained according to a back−transmission error algorithm, and Figure 3 schematically shows the structure of the neural network of the input layer, hidden layer, and output layer. Each layer contains many neurons inside, and the computational unit in the neuron is the activation function.  The network contains two stages in the calculation, the phase of forwardl sion of the signal and backwardly transmission of the error. Before training network, the target value of the final output result should be set first. In the training the neural network, the size of the weights and thresholds are constan and the actual value of the output is constantly close to the preset target valu optimal neural network was finally obtained.
The first step is the web forward transmission. Let the nodes in the input network be numbered j (j = 1, 2, …, v), the nodes in the hidden layer of the n assigned numbers n (n = 1, 2, …, u), and the nodes in the output layer of the n given numbers k (k = 1, 2, …, z). The signal passes the data through the input hidden layer. The input Nn of the nth node of the input layer is given by Equat where bn is the threshold of the nth node in the hidden layer, wnj is the weigh node in the input layer to the nth node in the hidden layer, and xj is the input v jth node in the input layer. Substituting Equation (1) into the activation functi hidden layer, the output of the nth node of the hidden layer is given as follow Using the output of the hidden layer as the input of the output layer, the i kth node in the output layer is obtained as follows: where ak is the threshold value of the kth node of the output layer, wkn is the w nth node in the input layer to the kth node in the output layer, and yn is the ou of the nth node of the hidden layer. Finally, the output equation for the kth n The network contains two stages in the calculation, the phase of forwardly transmission of the signal and backwardly transmission of the error. Before training the neural network, the target value of the final output result should be set first. In the process of training the neural network, the size of the weights and thresholds are constantly revised, and the actual value of the output is constantly close to the preset target value, and the optimal neural network was finally obtained.
The first step is the web forward transmission. Let the nodes in the input layer of the network be numbered j (j = 1, 2, . . . , v), the nodes in the hidden layer of the network are assigned numbers n (n = 1, 2, . . . , u), and the nodes in the output layer of the network are given numbers k (k = 1, 2, . . . , z). The signal passes the data through the input layer to the hidden layer. The input N n of the nth node of the input layer is given by Equation (1): where b n is the threshold of the nth node in the hidden layer, w nj is the weight of the jth node in the input layer to the nth node in the hidden layer, and x j is the input value of the jth node in the input layer. Substituting Equation (1) into the activation function f h in the hidden layer, the output of the nth node of the hidden layer is given as follows: Using the output of the hidden layer as the input of the output layer, the input of the kth node in the output layer is obtained as follows: where a k is the threshold value of the kth node of the output layer, w kn is the weight of the nth node in the input layer to the kth node in the output layer, and y n is the output value of the nth node of the hidden layer. Finally, the output equation for the kth node of the output layer was obtained by substituting the Equation (3) into the activation function f o of the output layer: The second step is the network back transmission. When the forward transmission of the signal was completed, the network will calculate the output error of the neurons in each layer in reverse from the output layer. Then the error gradient descent method was used to adjust the weights and thresholds of each layer. Firstly, the Mean Square Error (MSE) E p , of a single sample datum p was calculated by the following equation: where T k is the expected value of the kth node. Using the gradient descent method to find the weight and threshold correction amount of the output layer and the hidden layer, respectively. The specific equations are shown in Equations (6)-(9), where ∆w kn and ∆w nj represent the weight correction of the output layer and the hidden layer, respectively, ∆a k and ∆b n denote the threshold correction of the output layer and the hidden layer, respectively.
where η is the learning rate of the network, and P is the number of training samples.

Genetic Algorithm
GA is to generate a new population by continuously searching to select the optimal solution, the specific process is similar to the process of human genetic evolution. The specific process of GA is the selection and initialization of the population, fitness function, selection, crossover, mutation, getting new individuals, and judging whether the new individuals meet the required conditions. If they do not meet the conditions, then return to the selection operator to continue the search for the optimal solution. If they meet the conditions, then the optimal solution is the output. Individuals are often encoded by binary encoding, multiparameter cascade encoding, and floating−point encoding. The encoding methods have their characteristics and were adapted to different practical application scenarios. During the operation process, the individuals were selected and crossed between two operators. The individuals themselves are affected by the variation operators.

Computational Procedure of Back Propagation Neural Network Algorithm Based on Genetic Algorithm Optimization
In the present study, the process of the BP neural network algorithm based on GA optimization has four steps. Figure 4 shows the flow chart of the computation. In the first step, a dataset was created, and each joint angle was divided equally into 180 parts in its range of motion. The corresponding end location matrices were solved using the forward kinematic equations, and these matrices were transformed to correspond joint matrices. The 180 datasets were generated and stored as the elements of the input and output layers in the neural network, respectively. In the second step, a three−layer BP neural network model was built. The number of input layer neurons I was set to 12, the number of output layer neurons O was set to 7, and the number of hidden layer neurons H was calculated according to the empirical Equation (10). The number of neurons in the hidden layer was calculated to be 14 and the Pureline function [28] was selected as the activation function for the input and output layers in this study. The training times of the network were set to 1000, the training target was set to 0.0001, and the learning rate was set to 0.1.
where α is a regulation constant with a default value range of 1 to 10. The third step is to establish a GA model. This model used binary coding, the range of the weight threshold value was limited between −5 to 5, and the individual precision was set to 0.01. According to Equation (11), L was 9.967, and was rounded to 10, thus the individual binary coding length was 10. The primitive population was then determined, and the number of weights was equal to the number of neurons in the input layer multiplied by the number of neurons in the hidden layer, plus the number of neurons in the hidden layer multiplied by the number of neurons in the output layer. The number of thresholds was equal to the number of neurons in the hidden layer, plus the number of neurons in the output layer.
The network contains 287 weights and thresholds, of which 266 were weights and 21 were thresholds on the basis of the number of neurons in the three layers of the neural network. Forty groups of individuals were randomly generated as the parents and a matrix of size 40 × 2870 was formed and then was substituted into the neural network. The errors calculated by the neural network after bringing in the matrix were sorted to create the fitness function. Finally, the individuals were updated iteratively by three operators of selection, crossover, and variation until the set number of genetic generations was completed. As seen in Figure 5, with continuous evolutionary inheritance, the two norm error decreases In the first step, a dataset was created, and each joint angle was divided equally into 180 parts in its range of motion. The corresponding end location matrices were solved using the forward kinematic equations, and these matrices were transformed to correspond joint matrices. The 180 datasets were generated and stored as the elements of the input and output layers in the neural network, respectively. In the second step, a three−layer BP neural network model was built. The number of input layer neurons I was set to 12, the number of output layer neurons O was set to 7, and the number of hidden layer neurons H was calculated according to the empirical Equation (10). The number of neurons in the hidden layer was calculated to be 14 and the Pureline function [28] was selected as the activation function for the input and output layers in this study. The training times of the network were set to 1000, the training target was set to 0.0001, and the learning rate was set to 0.1.
where α is a regulation constant with a default value range of 1 to 10. The third step is to establish a GA model. This model used binary coding, the range of the weight threshold value was limited between −5 to 5, and the individual precision was set to 0.01. According to Equation (11), L was 9.967, and was rounded to 10, thus the individual binary coding length was 10. The primitive population was then determined, and the number of weights was equal to the number of neurons in the input layer multiplied by the number of neurons in the hidden layer, plus the number of neurons in the hidden layer multiplied by the number of neurons in the output layer. The number of thresholds was equal to the number of neurons in the hidden layer, plus the number of neurons in the output layer. The network contains 287 weights and thresholds, of which 266 were weights and 21 were thresholds on the basis of the number of neurons in the three layers of the neural network. Forty groups of individuals were randomly generated as the parents and a matrix of size 40 × 2870 was formed and then was substituted into the neural network. The errors calculated by the neural network after bringing in the matrix were sorted to create the fitness function. Finally, the individuals were updated iteratively by three operators of selection, crossover, and variation until the set number of genetic generations was completed. As seen in Figure 5, with continuous evolutionary inheritance, the two norm error decreases in steps, and finally, the weights and thresholds with the smallest errors are output.
where L is the number of binary codes, a is the lower bound for individual values, b is the upper bound on the value taken by an individual, and esp is the accuracy of the individual. In the fourth step, the optimal weights and thresholds obtained by the GA w stituted into the built BP neural network model. The first 150 sets of the dataset we for the training of the network, and the last 30 sets were used to verify the netwo ing results. As seen in Figure 6, with the number of training sessions increasing, t decreased continuously. The curve of the train and test set has almost the same tre the network meets the requirement. To verify the error of the trained network, the angles of the seven joints were into 500 parts and substituted into the trained network to analyze the error. The ma error between the solved data and the original data is 0.002 rad. As shown in Tabl 2, 10 sets of angle information were selected in the table. One was the original rad taset, and the other was the radian dataset solved by the trained network. By com the original radian data and the radian data, the error was very small. After verif training results of the network, the collected data of eight critical points were inpu trained network. The network outputs the angles of each joint, thus completing th lation of the inverse kinematic solution of the 7DOF robotic arm In the fourth step, the optimal weights and thresholds obtained by the GA were substituted into the built BP neural network model. The first 150 sets of the dataset were used for the training of the network, and the last 30 sets were used to verify the network training results. As seen in Figure 6, with the number of training sessions increasing, the MSE decreased continuously. The curve of the train and test set has almost the same trend, thus the network meets the requirement. In the fourth step, the optimal weights and thresholds obtained by the GA were sub stituted into the built BP neural network model. The first 150 sets of the dataset were used for the training of the network, and the last 30 sets were used to verify the network train ing results. As seen in Figure 6, with the number of training sessions increasing, the MSE decreased continuously. The curve of the train and test set has almost the same trend, thu the network meets the requirement. To verify the error of the trained network, the angles of the seven joints were divided into 500 parts and substituted into the trained network to analyze the error. The maximum error between the solved data and the original data is 0.002 rad. As shown in Tables 1 and 2, 10 sets of angle information were selected in the table. One was the original radian da taset, and the other was the radian dataset solved by the trained network. By comparing the original radian data and the radian data, the error was very small. After verifying th training results of the network, the collected data of eight critical points were input to the trained network. The network outputs the angles of each joint, thus completing the calcu lation of the inverse kinematic solution of the 7DOF robotic arm  To verify the error of the trained network, the angles of the seven joints were divided into 500 parts and substituted into the trained network to analyze the error. The maximum error between the solved data and the original data is 0.002 rad. As shown in Tables 1 and 2, 10 sets of angle information were selected in the table. One was the original radian dataset, and the other was the radian dataset solved by the trained network. By comparing the original radian data and the radian data, the error was very small. After verifying the training results of the network, the collected data of eight critical points were input to the trained network. The network outputs the angles of each joint, thus completing the calculation of the inverse kinematic solution of the 7DOF robotic arm We found the following methods for solving the inverse kinematic solution of the robotic arm in other articles: Jacobian transpose [29], neural network [30], Multiple Population Genetic Algorithm (MPGA) [31], Firefly Algorithm (FA), and Artificial Bee Colonies (ABC) [32]. The results of the comparison between these methods and the methods in this study are shown in Table 3. The comparison results in the table showed that the method used in this study had a higher accuracy than the other methods of calculation. Table 3. Comparison of robotic arm inverse kinematics solution methods.

Trajectory Planning of Redundant Robotic Arm for Upper Limb Rehabilitation
In this section, the rehabilitation curve was planned by the cubic B−spline interpolation method. The B−spline curve has the qualities of local support, composite and convex envelope compared with other curves, which can quickly obtain the calculation results. Thus, the cubic B−spline interpolation method was used to plan the rehabilitation trajectory. The inverse kinematic solution of the eight critical points was solved by the above GA−optimized BP neural network algorithm. The critical points were planned into a complete rehabilitation curve by B−spline interpolation in the Cartesian coordinate system, and the kinematic analysis of each joint was performed in the joint space to verify the smoothness of the trajectory motion.

Derivation of Cubic B−Spline Interpolation
During trajectory planning using the cubic B−spline interpolation method, the set of position points passed by each joint was known, and the planning of robotic arm trajectory was realized by back−calculating the control vertices of the B−spline curve to find the angle variables of each joint as a function of time.
The equation of the cubic B−spline curve function for each joint of the robot arm with the angle variable θ at time t is shown as follows: The first−order and second−order derivatives of Equation (12) could be obtained for the joint angular speed and joint angular acceleration for the time in Equation (13), where the value of time t ranges from 0 to 1, A is the basis function, V is the control vertex. The variable i (i = 1, 2, . . . , m) indicated the ith section of the curve, whose corresponding control vertices were numbered V i−1 , V i , V i+1 , V i+2 , and four adjacent control vertices constituted a set of control points.
From the previous sections, the B−spline curve has continuity, and the end of the ith segment of the curve was continuous with the start of the (i + 1)th segment of the curve, which yields Equation (14): Substituting Equation (14) into Equation (12), Equation (15) was obtained: From the continuity of the first−order and second−order derivatives of the B−spline curve, the following equation could be derived: Substituting Equation (16) into Equation (13) yielded the following results: A 0 (0), . .
Furthermore, according to the normality of the B−spline curve, the following equation could be obtained: In this study, the cubic B−spline curve interpolation method was used, so the equation basis function was set as follows, where a i , b i , c i , and d i are the four coefficients assumed: Combining Equations (12) to (19), an expression for the cubic B−spline curve equation for the joint angle at a given time in the joint space of the robot arm was obtained: (20) By solving the first and second−order derivatives of Equation (20), the equations of angular velocity and angular acceleration was obtained as follows: ..
To clarify the relation between the variables of Equation (20), it was rewritten in the form of a matrix with the following expression: The value of the control vertex V i was unknown, and the control vertex was found by inputting the joint coordinates of the robotic arm. The points on the trajectory of the actuators were known, and the joint angle was obtained by the inverse kinematics. Taking joint #1 as an example, the curve passed through points P 1 , P 2 , . . . , P m , and the two adjacent points were planned by a section of the cubic B−spline curve. A total of m − 1 segments of the B−spline curve existed. In addition, because of the continuity of the curve, Equation (24) could be obtained: The equation containing m data points described the geometry of the curve or surface and m + 2 control point vertices. To ensure the smoothness and continuity of the curve, two constraints were added as V 0 = V m and V 1 = V m+1 . Thus Equation (24) was transformed into the form of a matrix: The system of equations had m equations and unknown quantities, which was solved for a unique set of solutions. The values of the control vertices V 0 to V m−1 was obtained by this solution set. Subsequently, all the control vertices V 0 to V m+1 were derived by the two constraints. Thus, the solution was completed for each section of the trajectory B−spline curve equation.

Rehabilitation Trajectory Planning under Cartesian Space
Under the Cartesian space, the coordinates of the eight critical points, which were solved by the inverse kinematic solution in Section 2, were transformed to the coordinates under the Cartesian coordinate system by the forward kinematic solution. The trajectories of the eight critical points after the forward kinematic solution were planned using the cubic B−spline interpolation method derived above to form a complete section of the spatial curve, as shown in Figure 7. The green line in the figure is the original rehabilitation trajectory. The blue line is the optimized trajectory with an increasing movement range. For people with upper limb disorders, there are many unreachable spatial locations. In this study, we assumed that the patient had a disability in the left upper limb, and the left side of the patient corresponds to the negative direction of the Y−axis in the figure. However, after a period of upper limb rehabilitation, the patient's reachable range of the upper limb increased. Therefore, the motion trajectory of the robotic arm needs to expand accordingly, which is the reason why the trajectory optimization was needed.
nsors 2022, 22, x FOR PEER REVIEW limb increased. Therefore, the motion trajectory of the robotic arm nee cordingly, which is the reason why the trajectory optimization was need

Joint Kinematic Analysis of Redundant Robotic Arm for Upper Limb Rehabilitation
Based on the trajectory planned with the cubic B−spline interpolatio the Cartesian space above, the kinematic analysis was performed for eac joint space to verify the feasibility of the trajectory. Take joint #1 as an ex at each of the pre−given critical points during the motion were: q = [81.40° 117.92°, 82.55°, 42.88°, 38.98°, 42.42°]. Set the corresponding time interval angle difference between adjacent points, specific time intervals were: t s, 4.0 s, 4.0 s, 2.0 s, 2.0 s] and it took a total of 25 s to complete the rehabili As shown in Figure 8, the circles in the location curve represent the eig and the cubic B−spline interpolation method is used to connect the eigh We wrote an algorithmic program based on MATLAB to calculate the sp ation of the joints. For the acceleration calculation, we only calculated th several locations during the motion and connected them with straight lin we mainly focus on the continuity of the speed curve and the curves pas points and the cubic B−spline interpolation method was able to adapt to t motion and completed the planning of the rehabilitation trajectory. The tained smoothly during the motion, which verified the stability of the r botic arm motion.

Joint Kinematic Analysis of Redundant Robotic Arm for Upper Limb Rehabilitation
Based on the trajectory planned with the cubic B−spline interpolation method under the Cartesian space above, the kinematic analysis was performed for each joint under the joint space to verify the feasibility of the trajectory. Take joint #1 as an example, its angles at each of the pre−given critical points during the motion were:  Figure 8, the circles in the location curve represent the eight critical points, and the cubic B−spline interpolation method is used to connect the eight critical points. We wrote an algorithmic program based on MATLAB to calculate the speed and acceleration of the joints. For the acceleration calculation, we only calculated the acceleration at several locations during the motion and connected them with straight lines. In this study, we mainly focus on the continuity of the speed curve and the curves passed all the preset points and the cubic B−spline interpolation method was able to adapt to the complex joint motion and completed the planning of the rehabilitation trajectory. The speed was maintained smoothly during the motion, which verified the stability of the rehabilitation robotic arm motion. In addition to interpolation method of the cubic B−spline curve, we also used th cubic polynomial interpolation, the five−polynomial interpolation, and the seven−polyno mial interpolation method to plan the motion trajectory of the robotic arm. In addition we performed trajectory planning as well as MATLAB simulation for each method fol lowing the procedure of cubic B−spline interpolation.
First is the cubic polynomial interpolation method, the simulation results are show in Figure 9. The second is the five−polynomial interpolation method, the simulation results ar shown in Figure 10. In addition to interpolation method of the cubic B−spline curve, we also used the cubic polynomial interpolation, the five−polynomial interpolation, and the seven−polynomial interpolation method to plan the motion trajectory of the robotic arm. In addition, we performed trajectory planning as well as MATLAB simulation for each method following the procedure of cubic B−spline interpolation.
First is the cubic polynomial interpolation method, the simulation results are shown in Figure 9. In addition to interpolation method of the cubic B−spline curve, we also used the cubic polynomial interpolation, the five−polynomial interpolation, and the seven−polyno mial interpolation method to plan the motion trajectory of the robotic arm. In addition we performed trajectory planning as well as MATLAB simulation for each method fol lowing the procedure of cubic B−spline interpolation.
First is the cubic polynomial interpolation method, the simulation results are shown in Figure 9. The second is the five−polynomial interpolation method, the simulation results are shown in Figure 10. The second is the five−polynomial interpolation method, the simulation results are shown in Figure 10. Finally, the simulation results of the seven−polynomial interpolation are shown in Figure 11. The simulation results show that the location and speed of the trajectories planned by these three methods were relatively stable and ensured the smoothness of the motion of the robotic arm. Finally, the simulation results of the seven−polynomial interpolation are shown in Figure 11. Finally, the simulation results of the seven−polynomial interpolation are shown i Figure 11. The simulation results show that the location and speed of the trajectories planne by these three methods were relatively stable and ensured the smoothness of the motio of the robotic arm. The simulation results show that the location and speed of the trajectories planned by these three methods were relatively stable and ensured the smoothness of the motion of the robotic arm.

Optimization of Rehabilitation Trajectories with the Goal of Time Minimization
Since different patients have different needs for the upper limb rehabilitation, we need to expand the applicable range of the redundant robotic arm for upper limb rehabilitation. We can change certain parameters in the rehabilitation process using trajectory optimization to achieve the purpose of expanding the applicable range. For example, we could change the motion time of the end of the robotic arm between the critical points and thus change the speed of the end movement. We could also modify the location of the critical points and change the shape of the rehabilitation trajectory to meet the needs of different patients for rehabilitation. The curve after trajectory planning consists of multiple cubic B−spline curves that have been connected. Taking the total running time of the robot arm as the optimization objective, the running time of each segment of the curve was set to t i . The time of each part of the trajectory was optimized separately to obtain the final result. The objective function was set as follows: where T is the total operating time of the robotic arm. For time−targeted optimization, the motion time could not be reduced indefinitely and constraints needed to be added. Kinematic constraints include speed, acceleration, and jerk (derivative of acceleration). First is the speed constraint. The joint angular speed equation was derived from the joint equation to the first−order derivative of time, the maximum angular speed of the robotic arm in motion might be at any point on the interval when the robotic arm motion time was at t i to t i+1 , there was an expression as follows: where . θ is the maximum angular speed, . θ i is the absolute value of the angular speed at moment t i , and . θ ix is the maximum angular speed in the interval t i to t i+1 . It was difficult to solve the equation for the absolute maximum value, and the golden mean method was used to solve the equation for the maximum value of . θ max . The constraint form was given as follows: The second is the acceleration constraint. The joint angular acceleration equation was obtained from the joint angle equation by taking the second−order derivative of time and solving for the maximum value in the same way as the speed constraint. The maximum joint angular acceleration equation was solved as follows: ..
θ is the maximum angular acceleration, ..
θ i is the absolute value of angular acceleration at moment t i , and .. θ ix is maximum angular acceleration in absolute value in the interval t i to t i+1 . The acceleration constraint form of the robotic arm was as follows: The third is the jerk constraint. The joint angular jerk equation was solved for the maximum joint angular acceleration equation by taking the third−order derivative of the joint angle equation for a time and was given as follows: where ... θ is the maximum angular jerk, | ... θ i | is add the absolute value of angular jerk at moment t i , and | ... θ ix | is the absolute maximum plus angular jerk in the interval t i to t i+1 .

MATLAB Numerical Simulation of Redundant Robotic Arm for Upper Limb Rehabilitation
We assumed that the model had m joint angle values and the motion time between each adjacent two points was t i , a GA was used to find the optimal solution for t i in the model. In this study, the time interval matrix between two adjacent points was set as t = [8.0 s, 2.5 s, 2.5 s, 4.0 s, 4.0 s, 2.0 s, 2.0 s, and 50 groups of individuals were randomly generated at each given time. For each random individual, the fitness size f was assigned by the following equation.
When the constraint was satisfied The final results were calculated by the iterating depending on the iteration rules of the GA. According to the application scenario of the upper limb rehabilitation, the maximum joint angular speed V max in the algorithm was set to 15 • /s, the maximum joint angular acceleration A max was set to 20 • /s 2 , and the maximum joint angular jerk I max was set to 100 • /s 3 . Under the condition of satisfying the kinematic constraints, the total motion time of the robotic arm before the optimization was 25 s, which was reduced to 20.44 s after the time optimization. The time comparison before and after track optimization is shown in Table 4. We wrote the above algorithm based on MATLAB. The kinematic analysis was performed on the curve obtained by cubic B−spline interpolation for joint #1 of the robotic arm and the trajectory optimized by the GA. As shown in Figure 12, the red line indicates the unoptimized analysis results and the blue line represents the optimized results. The results show that after optimization, the location and angular speed profiles of the joints were smooth, and the rehabilitation movements were completed in a shorter period of time, which expanded the range of applicability. In terms of safety, the patient had a certain tolerance after a period of rehabilitation training. Only after the patient had received a period of formal rehabilitation training, the time was shortened and the rehabilitation effect was improved. In addition, even though the acceleration curve changed more than the unoptimized curve, this change was still in a controlled range to ensure that the patient can receive secondary injuries during the rehabilitation process. We have read the literature to the extent of our ability and describe the structure and movement of the upper limb as follows [33].
In terms of the upper limb structure, the upper limb consists of the shoulder girdle and the arm. The shoulder girdle consists of the clavicle and the scapula. The shoulder joint, which connects the arm to the scapula, is a multiaxial synovial ball and socket joint. The humeral head was shaped as a ball which coincides with the spherical surface of the glenoid cavity of the scapula, and the surfaces of the humeral head and cavity were covered with articular cartilage. In addition, the joint capsule was strengthened by muscles, tendons, and ligaments such as the glenohumeral, coracohumeral, and coracoacromial ligaments.
The arm was built of the humerus, ulna, radius bone, 8 carpal bones of the wrist, 5 metacarpal bones, and 14 phalanges of the fingers. The elbow consists of the humeroulnar, humeroradial, and proximal radioulnar joints, which were surrounded by a common joint capsule. The frontal stability was mainly ensured by the ulnar collateral ligament and the radial collateral ligament. Patients with impaired upper limb function usually suffer from wear and tear of the ligaments or joints of the upper limb. Therefore, the upper limb rehabilitation robotic arm in this study treated the patient's upper limb by helping the patient to perform upper limb movements through end−traction method of the robotic arm.
In terms of the upper limb movement, the upper limb skeletal movement was produced by the 43 muscles being connected to structures mentioned above. The shoulder joint allows for abduction/adduction, flexion/extension, and internal/external rotation. At the elbow joint, there is elbow flexion/extension and forearm supination/pronation. In the wrist, there is flexion/extension and radial/ulnar deviation. However, because of the presence of the upper limb joints, there is a limited range of motion in each of these joints. Therefore, the upper limb rehabilitation of the robotic arm in this study planned the rehabilitation trajectory by pre−given critical points, thus avoiding the rehabilitation trajectory to reach the position where the patient's upper limb cannot reach, which can ensure the safety of the patient during the rehabilitation process.
This redundant robotic arm for upper limb rehabilitation is still in the stage of theoretical modeling, numerical simulation, development, and testing of the prototype. Therefore, for safety reasons, it has not been commissioned and applied to patients yet. However, to obtain a better perceptual feedback on the human body, we had conducted preliminary tests with the manufactured robotic arm on several patients with healthy upper limb. These patients, with their sensitive upper limb perception, were able to provide better feedback on the treatment effect of the upper limb rehabilitation of the robotic arm, the We have read the literature to the extent of our ability and describe the structure and movement of the upper limb as follows [33].
In terms of the upper limb structure, the upper limb consists of the shoulder girdle and the arm. The shoulder girdle consists of the clavicle and the scapula. The shoulder joint, which connects the arm to the scapula, is a multiaxial synovial ball and socket joint. The humeral head was shaped as a ball which coincides with the spherical surface of the glenoid cavity of the scapula, and the surfaces of the humeral head and cavity were covered with articular cartilage. In addition, the joint capsule was strengthened by muscles, tendons, and ligaments such as the glenohumeral, coracohumeral, and coracoacromial ligaments.
The arm was built of the humerus, ulna, radius bone, 8 carpal bones of the wrist, 5 metacarpal bones, and 14 phalanges of the fingers. The elbow consists of the humeroulnar, humeroradial, and proximal radioulnar joints, which were surrounded by a common joint capsule. The frontal stability was mainly ensured by the ulnar collateral ligament and the radial collateral ligament. Patients with impaired upper limb function usually suffer from wear and tear of the ligaments or joints of the upper limb. Therefore, the upper limb rehabilitation robotic arm in this study treated the patient's upper limb by helping the patient to perform upper limb movements through end−traction method of the robotic arm.
In terms of the upper limb movement, the upper limb skeletal movement was produced by the 43 muscles being connected to structures mentioned above. The shoulder joint allows for abduction/adduction, flexion/extension, and internal/external rotation. At the elbow joint, there is elbow flexion/extension and forearm supination/pronation. In the wrist, there is flexion/extension and radial/ulnar deviation. However, because of the presence of the upper limb joints, there is a limited range of motion in each of these joints. Therefore, the upper limb rehabilitation of the robotic arm in this study planned the rehabilitation trajectory by pre−given critical points, thus avoiding the rehabilitation trajectory to reach the position where the patient's upper limb cannot reach, which can ensure the safety of the patient during the rehabilitation process.
This redundant robotic arm for upper limb rehabilitation is still in the stage of theoretical modeling, numerical simulation, development, and testing of the prototype. Therefore, for safety reasons, it has not been commissioned and applied to patients yet. However, to obtain a better perceptual feedback on the human body, we had conducted preliminary tests with the manufactured robotic arm on several patients with healthy upper limb. These patients, with their sensitive upper limb perception, were able to provide better feedback on the treatment effect of the upper limb rehabilitation of the robotic arm, the comfort level during treatment (no sudden changes in force) and the continuity of the rehabilitation trajectory (continuous and smooth location and speed of the robotic arm). The rehabilitation movement of the end traction was performed according to the planned trajectory in this study, as shown in Figure 13. rehabilitation trajectory (continuous and smooth location and speed of the robotic arm). The rehabilitation movement of the end traction was performed according to the planned trajectory in this study, as shown in Figure 13. The rehabilitation trajectory passed the eight critical points mentioned in this study, and the positions and postures of the robotic arm and the treated person during the rehabilitation process corresponded to a to h in Figure 13. The feedback from the treated person could better reflect the rehabilitation effect of the robotic arm, so that the trajectory of the robotic arm could be adjusted and optimized accordingly before practical application. For example, the trajectory optimization method with the goal of time minimization in this study could change the speed during the motion and the magnitude of the force applied to the upper limb. Alternatively, the range and accuracy of the motion trajectory of the robotic arm end could be adjusted by changing the location of the critical points and adjusting the number of critical points. The test results show that the robotic arm did not cause excessive compression on the upper limb during the rehabilitation process guided by the end. The force and speed applied to the human upper limb were moderate, and the end of robotic arm was relatively smooth during the movement, which could achieve the effect of assisting in guiding the upper limb rehabilitation.
We will cooperate with the TCM rehabilitators from Beijing Red Medical Star Intelligent Technology Development Company Limited in the future and use the robotic arm in the actual rehabilitation process. We will discuss the effectiveness of the treatment and the practical application in the next articles.

Conclusions
In this study, we proposed a BP neural network algorithm optimized by GA. Compared with other algorithms, the computational results of this algorithm have better accuracy. We used cubic B−spline interpolation method for trajectory planning and per- The rehabilitation trajectory passed the eight critical points mentioned in this study, and the positions and postures of the robotic arm and the treated person during the rehabilitation process corresponded to a to h in Figure 13. The feedback from the treated person could better reflect the rehabilitation effect of the robotic arm, so that the trajectory of the robotic arm could be adjusted and optimized accordingly before practical application. For example, the trajectory optimization method with the goal of time minimization in this study could change the speed during the motion and the magnitude of the force applied to the upper limb. Alternatively, the range and accuracy of the motion trajectory of the robotic arm end could be adjusted by changing the location of the critical points and adjusting the number of critical points. The test results show that the robotic arm did not cause excessive compression on the upper limb during the rehabilitation process guided by the end. The force and speed applied to the human upper limb were moderate, and the end of robotic arm was relatively smooth during the movement, which could achieve the effect of assisting in guiding the upper limb rehabilitation.
We will cooperate with the TCM rehabilitators from Beijing Red Medical Star Intelligent Technology Development Company Limited in the future and use the robotic arm in the actual rehabilitation process. We will discuss the effectiveness of the treatment and the practical application in the next articles.

Conclusions
In this study, we proposed a BP neural network algorithm optimized by GA. Compared with other algorithms, the computational results of this algorithm have better accuracy. We used cubic B−spline interpolation method for trajectory planning and performed numerical simulation with MATLAB. The results show that the location and speed profiles of the planned curve were relatively smooth and could meet the needs of patient rehabilitation. To further improve the rehabilitation effect, we optimized the planned trajectory based on GA with the goal of time minimization, and the numerical simulation results in MATLAB showed that the optimized curve reduced the motion time, thus increasing the motion speed of the end of the robotic arm and improving the rehabilitation effect. For safety reasons, we have not yet applied this upper limb rehabilitation of the robotic arm to patient rehabilitation treatment. We had tested it on several patients with normal upper limb and the results show that the force and speed applied to the upper limbs during the rehabilitation process was relatively stable and could achieve the effect of assisting in guiding the rehabilitation of the upper limbs. In the future, we will cooperate with TCM rehabilitators and use the upper limb rehabilitation of the robotic arm in the rehabilitation process of patients. We will discuss the treatment effect and practical application of the upper limb rehabilitation arm in our future work.